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Abstract. We compute the spectral density for ensembles of of sparse symmetric 
random matrices using replica, managing to circumvent difficulties that have been 
encountered in earlier approaches along the lines first suggested in a seminal paper by 
Rodgers and Bray. Due attention is payed to the issue of localization. Our approach is not 
restricted to matrices defined on graphs with Poissonian degree distribution. Matrices 
defined on regular random graphs or on scale-free graphs, are easily handled. We also 
look at matrices with row constraints such as discrete graph Laplacians. Our approach 
naturally allows to unfold the total density of states into contributions coming from 
vertices of different local coordination. 



1. Introduction 

Since its inception by Wigner in the context of describing spectra of excited nuclei [1], 
Random Matrix Theory (RMT) has found applications in numerous areas of science, 
including questions concerning the stability of complex systems [2j, electron localisation 
[3], quantum chaos [4], Quantum Chromo Dynamics [5], finance [6j [7j, the physics of 
glasses both at elevated [8j [9] and low [TUl E] temperatures, number theory [12J, and 
many many more. For an extensive review describing many of the applications in physics 
see, e.g. [13]. 

In the present paper we revisit the problem of determining the spectral density for 
ensembles of sparse random matrices pioneered two decades ago in seminal papers by 
Bray and Rodgers [T4j lT5"]. The problem has in recent years received much renewed interest 
in connection with the study of complex networks, motivated, for instance, by the fact 
that geometric and topological properties of networks are reflected in spectral properties 
of adjacency matrices defining the networks in question JT6J [T7j. Also, phenomena such 
as non-exponential relaxation in glassy systems and gels [151 EE] - intimately related 
to Lifshitz tails [19] and Griffiths' singularities in disordered systems [20] - as well as 
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Anderson localization of electronic [21] or vibrational [22] states have been studied in 
sparsely connected random systems, as finite dimensional versions of these problems have 
proven to be extremely difficult to analyse. A wealth of analytical and numerical results 
has been accumulated on these systems in recent years. Progress has, however, been 
partly hampered by the fact that full solutions of the Rodgers-Bray integral equation 
|14j . in terms of which spectral densities of the sparse random matrices in question are 
computed, have so far eluded us. Asymptotic analyses for large average connectivities 
[H1[T5], an d other approximation schemes such as the single defect approximation (SDA) 
and the effective medium approximation (EMA) [231 I2H EH] or very recently [25], as well 
as numerical diagonalization (e.g. [26]) had to come in for help. 

In what follows we describe some significant progress in the understanding of this problem, 
based upon advances in the statistical mechanical analysis of sparsely connected spin- 
glass like systems seen in the last couple of years [23 [2H] - in the present context in 
particular the proposal of a stochastic population-dynamics algorithm [28j to solve the 
nonlinear integral equations appearing in the solution of these problems, and the recent 
generalization of these methods to systems with continuous degrees of freedom, such as 
models of sparsely connected vector spins [29], or finitely coordinated models for low- 
temperature phases of amorphous systems [30] . 

It is well known that the average spectral density of an ensemble M. of iV x iV matrices 
M can be computed from the ensemble average of the imaginary part of their resolvent 
via 

JnW = ^Im Tr [A.I-M]- 1 , (1) 

in which II is the N x N unit matrix, and A e = A — ie, the limit e — > + being understood. 
Following Edwards and Jones [31] , one can express this result in terms of the Gaussian 
integral 

f N dui f i } 

Z N = / II r—T 6XP 1 ~ o ^ u ^ X ^j ~ M a) u 3 ) ( 2 ) 

J i=l \jlTlll { * i,j J 

as 

2 d 1 N 

Piv(A) = -— Im — In Z N = -Re £ («?) , (3) 

y '' i=l 

using the replica method to evaluate the average of the logarithm in ([3]) over the ensemble 
M. of matrices M under consideration. The 'averages' (uf) in ([3]) are evaluated with 
respect to the 'Gaussian measure' defined by ([2])Jj3 This has been the path taken in [T3] ; 
we shall initially follow their reasoning. 

| Note that we are using probabilistic notions in a loose, metaphorical sense, as the Gaussian measures 
used in these calculations are complex. 
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Disregarding the complex nature of the 'Hamiltonian' in the evaluation of (121) , the 
mathematical problem posed in (J2J), ([3]) is analogous to the evaluation of an 'internal 
energy of a disordered system with quenched disorder. Within the general class of finitely 
coordinated amorphous model systems considered in [30], the one represented by (j2J), (EJ) 
constitutes a particular sub-class, viz. that of harmonically coupled systems, for which the 
analysis was found to be much simpler than for systems involving anharmonic couplings. 
Indeed, while the solution of the latter required the self-consistent determination of 
probability distributions over infinite dimensional function-spaces, it was realized in 
[30] that solutions of harmonically coupled systems could be formulated in terms of 
superpositions of Gaussians, and that the self-consistency problem reduced to the (much 
simpler) problem of a self-consistent determination of the probability distribution of their 
variances. 

It can be fairly argued that this last insight is, in fact, easier to obtain within a Bethe- 
Peierls or cavity type approach [28], in which (jSJ is recursively evaluated for given 
instances on graphs which are locally tree-like, ignoring correlations among subtrees 
- an approximation that becomes exact, e.g., for random graphs that remain finitely 
coordinated in the thermodynamic limit. This approach is taken in a separate publication 
[32] . in which (finite) single-instances and promising algorithmic aspects of the problem 
are being highlighted. 

Although [30] describes all technical details needed for a replica analysis of the present 
problem, we shall nevertheless reproduce the key steps here, both to keep the paper self- 
contained, and to point out along the way where the impasse in [T3] arises, and how it is 
circumvented. 

The remainder of the paper is organized as follows. In Sec. 2, we describe the replica 
analysis of the problem posed by (J2J) , (ED , specializing to matrices defined on Poissonian 
(Erdos-Renyi) random graphs. It has been known for some time [3T| IT4] that the replica- 
symmetric high-temperature solution — i.e., a solution preserving both, permutation- 
symmetry among replica, and rotational symmetry in the space of replica — is exact 
for problems of the type considered here. Accordingly, a representation that respects 
these symmetries is formulated in Sec. 2.1. It is at this point where our formulation 
departs from that of [14]. In Sec. 3 we present results for a variety of examples, 
and compare with numerical diagonalization results for large finite matrices to assess 
their quality. In sufficiently sparse graphs, one expects localized states to appear. The 
signatures of localization within our approach are discussed throughout Sec. 3, with 
inverse participation ratios (IPRs) as a diagnostic tool looked at in Sec. 3.2. A detailed 
investigation of Anderson localization for (discrete) Schrodinger operators on sparse 
random graphs will be reserved to a separate publication [33]. Matrices with bimodal 
instead of Gaussian random couplings are studied in Sec. 3.3. As the formal structure of 



Spectra of Sparse Random Matrices 



4 



the self-consistency problem remains unaltered when the Poissonian random graphs are 
replaced by graphs with other degree distributions [30] , we can exploit this fact to present 
results for regular and scale-free random graphs in Sec. 3.4. Modifications needed to treat 
matrices with row-constraints, such as discrete graph Laplacians are outlined in Sec. 3.5. 
Our approach naturally allows to unfold the total density of states into contributions 
coming from vertices of different local coordination, and we finally present an example of 
such an unfolding in Sec. 3.6. The final Sec. 4 contains a brief summary and an outlook 
on promising directions for future research. 



2. Replica Analysis 

2.1. General Formulation 

Here we briefly outline the evaluation of (|2J), (JSJ) for sparse symmetric matries M of the 
form 

My = <; J K, J , (4) 

in which C = {c{j} is a symmetric adjacency matrix of an undirected random graph 
(with ca = 0), and the non-zero elements of M are specified by the Kij, also taken to be 
symmetric in the indices. Within the present outline we restrict ourselves for the sake of 
simplicity to adjacency matrices of Erdos-Renyi random graphs, with 

p (i c ij}) = lI/'f'y)<V....... and p{cij) = ( 1 - -^j 5 Ctjfi + ^5 ClJ ,i , 



exhibiting a Posisssonian degree distribution with average coordination c. We note at the 
outset that formal results carry over without modification to other cases [30]. There is 
no need at this point to specify the distribution of the Ky, but we shall typically look at 
Gaussian and bimodal distributions. 



The average ([3]) is evaluated using replica \nZ N = lim„^ h hiZjy-, starting with integer 
numbers of replica as usual. After performing the average over the distribution of the 
connectivities one obtains 

^ = /n^^= ex p(-^ A ^I] M m + ^E (( ex P [iK^UiaU^ > ( 5 ) 



i,a ij \ \ \ o, /IK 



in which {■ ■ -)k refers to an average over the distribution of the K^. A decoupling of sites 
is achieved by introducing the replicated density 

p( u ) = jjJ2Ii s { Ua ~ Uia ) ' ( 6 ) 
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with u denoting the replica vector u = (ui,u 2 , ■ ■ ■ ,u n ), and enforcing its definition via 
functional 5 distributions, 

1 = J VpVp exp i-ij dup(u) (n P (u) - ElK^a - 11*)] I . (7) 

This gives (using shorthands of the form dp(u) = dup(u) where useful) 



Z n N = Jvpjvp exp jjV °- J dp(u)dp(v) (Jexp [iK^2u a v a j J 
- j du ip(u)p{u) + In J [] "4= exp f i p{x) - l - X £ J2 



- 1 



(8) 



allowing to evaluate N~ l \aZ^ by a saddle point method. The stationarity conditions 
w.r.t. variations of p and p read 



^)^/dp W ((e XP (,A-i:^: 



(9) 



K 



exp % p(u) - \ \ £ J2 a K 



p{u) = ^— ■ ^- . (10) 

Jduexp iip{u) -|A e E n «2 

The way in which sites are decoupled constitutes the first point of departure between our 
treatment and that of [H] and subsequent analyses inspired by it (e.g. [3^1 [35]). In these 
papers the averaged exponential expressions in the exponent of t^j, 

f(Ui ■ Vj) = f[Yl U i* V i°) = ( eX P ( iK Yl U iaVjc) ) - 1 , (H) 
a \ a Ik 

is expanded, and an infinite family of multi-replica generalizations of Edwards Anderson 
order parameters (and corresponding Hubbard- Stratonovich transformations) are used to 
decouple the sites, much as in the treatment of the dilute spin-glass problem by Viana 
and Bray [36] . The authors then use the expansion and the infinite set of self-consistency 
equations for the multi-replica generalizations of Edwards Anderson order parameters to 
construct a non-linear integral equation for a function g defined via a suitable 'average' 
of /; see [H] for details. Our treatment in this respect is closer in spirit to the alternative 
approach of Kanter and Sompolinsky [37] who treat local field distributions (which in the 
general context of disordered amorphous systems discussed in [30] become distributions 
of local potentials) as the primary object of their theory. 

However, the difference between our treatment and that of [2] is at this point still 
superficial. Indeed, we have the correspondence 



ip(u) = cg(u) 



(12) 
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between our 'conjugate density' p and the function g of [H|. With this identification, 
and ffTUj) can be combined to give 



which is the Rodgers-Bray integral equation for general distributions of non-zero bond 
strengths. 

2.2. Replica Symmetry 

To deal with the n — ► limit in these equations, assumptions concerning the invariance 
properties of the solutions p(u) and p(u) of (jUj) and ffTUl) — alternatively of the solution 
g(u) of (|T3|) — under transformations among the replica are required. It has been 
established for some time [3TL [H] that the replica-symmetric high-temperature solution 
- i.e., a solution preserving both, permutation-symmetry among replica, and rotational 
symmetry in the space of replica — is exact for problems of the type considered here. 
It is here where the paths taken in the present paper and in [H] really bifurcate. In 
|14j . the assumption g(u) = g(u), with u = \u\ is used to perform the angular integrals 
in n-dimensional polar coordinates in (TT3l) . resulting in an integral equation for g(u) in 
the n — > 0-limit. This integral equation has also been obtained using the supersymmetry 
approach [38] . It has, however, so far resisted exhaustive analysis or full numerical solution. 

In the present paper we follow [HD], and represent p and p as superpositions of replica- 
symmetric functions, using the observation made in [30] that superpositions of Gaussians 
of the form 



would provide exact solutions for harmonically coupled systems. Note that these 
expressions do indeed preserve permutation symmetry among replica as well as rotational 
symmetry. In (fT5|) the constant c is to be determined such that n is normalized, 
/ dfi(uj) = 1. We note that these representations make sense only for Re uj > and 
Re Cj > 0; later on we shall find that these conditions are self-consistently met for solutions 
of the fixed point equations. Expressing ([H]) in terms of 7r and it, we get 




(13) 




(15) 



(14) 



Z n N = J VTTVnex P {N[G 1 [TT]+G 2 [7r,n] + G 3 [n]]} 



(16) 
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As n — > 0, the functionals G\, G2 and G3 evaluate to 



Gi[/t] ~ n- I dir(uj)dii(uj') (in 



Z 2 (u,u',K) \ 
Z(u)Z(u') I 
Z(uj + uj) 



K 



G , 2 [7r,7r] ~ —c — ncj d7t(uj)dn(uj) In 
G 3 [tt] ^c + n /{d7T} fc ln-A 



Z(uj)Z(uj) 

Z x ( {uj} k ) 
Z(u t ) 



(17) 
(18) 
(19) 



in which we have introduced the shorthands {d7r}fc = Fu =1 d7r(a^), and {uj}k = Y$=i&h 
a Poissonian connectivity distribution 



Pc{k) 



k\ 



exp|— c| 



with average connectivity (k) — c, and the 'partition functions' 



Z(w) 

Z 2 (u,u',K) 



du exp 
dw 



— 

2 



^/27r/w , 



exp 



2n/i 
dudv exp 



~(«A e + {^} fc )u 2 



1/2 



2iKuv 



2tt 



(20) 

(21) 
(22) 

(23) 



y/uu' + K 2 

Note that the 0(1) contributions of G2 and G3 in the exponent of (jHJ) cancel in their sum. 

The stationarity condition of the functional integral (jSJ) w.r.t variations of p and p is 
reformulated in terms of stationarity conditions w.r.t variations 7r and 7r, 

Z(a) + uj) 



c J d.7r(a>) In 
dn(uj) In 



Z(u>)Z(u;) 
Z(a> + u) 



c I dn(u>) L Z ^'; K A +p 



Z uj )Z uj' 



K 



kpc(k) J {d7r} fe _i In 



Z\ e [uj + {uj}k-i) 
Z{u)Ht=}Z(u e ) 



(24) 



+ A(,25) 



with and /t Lagrange multipliers to take the normalization of 7r and 7r into account. 

The conditions that (1241) must hold for all uj and similarly that (125|) must hold for all uj 
can be translated |28| into 



tt(w) = - / d7r(u/) - ft(u/,.FQ) 



k>l C 



1T(UJ 



(26) 
(27) 



in which Cl(u',K) and f2({u>}fc_i) are defined via 

Z 2 (u,u',K) 



Z(uj + n(uj',K)) 



Z(u' 



K 2 



0J' 



(2* 
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and 

fc-i 



n({&} h - 1 ) = i\ e + Y,&t, (29) 



respectively. Given that n is normalized, it follows from ( 12611 that the same is true for ft, 
provided c = c, so the fixed point equations take their final form as 

tt(w)= Jdn(uj') (5(w-Cl(u/,K))) K , (30) 

7r(co) = J2 ~Pc(k) [ {dn} k ^ 5 (u - fl({w} w )) . (31) 
fc>i c 

These equations can be seen as special cases of the general framework derived in [30], 
when restricted to harmonically coupled random systems. In [30] it is shown that they 
hold — unmodified — for non-Poissonian degree distributions as well, as long as the 
average connectivity in these systems remains finite. 

Note that for all e > 0, ir and H — self-consistently - - have support in Re uj > 
and Re u > as required. The equations take a form that suggests solving them via a 
stochastic population-based algorithm, as described in Appendix A. 

For the thermodynamic limit of the spectral density we obtain from (j2J), (EJ) and ( fT6l) - (l23l) 
that 



1 oo 

p(A) = -Im 5> c (fc) /{dTT} 



i\ £ + {u} k 

= l -t Pc(k) /{dvr}, Re({ f fc + £) - 2 . (32) 

nfo J (Re(M fc + e)) 2 + (A + Im{cI ) } fc ) 2 1 ; 

This expression has a natural interpretation as a sum of contributions of local-densities of 

state of sites with connectivities k, weighted according to their probability of occurrence. 

Referring to (j3J), we may further identify the 

o\ = -Im — ^ (33) 

71 l\ £ + {uj} k 

as realizations of the variance of (Gaussian) marginals on sites of coordination k. 

With an eye towards disentangling singular (pure point) and continuous contributions to 
the spectral density, we find it useful to define 

P(a, b) = ^Pc(k) f {dn} k 5 (a - Re {u} k ) 5(b-lm {u} k ) , (34) 
k J 

with a > by construction. The density of states can then be expressed as an integral 
over P, 

da db , , . a + e 



, (A)=/ _ P(a , 6) ______ . (35) 
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Noting the singlular nature of the above integrand in the limit e — > for a = 0, we propose 
to isolate possible singular contributions to the spectral density by writing 

P(a,b) =P (b)6(a)+P(a,b) . (36) 

This gives 

/f (\n db ~ n A- f 

db P (b)£ E (b + A) + / P(a, b) V l(h+\V ' W 

J a >o 7T [a + e) z + (b + \) z 

in which C £ denotes a Lorentzian of width e. Our results below strongly suggest that, 
when the limit e — > is taken — thereby C £ (x) — > 5(x) — a non-zero value of 

P (-A) = lim / db P (b)C £ (b + A) (38) 

gives the contribution of the pure-point spectrum, originating from localized states, to 
the overall spectral density. 

This concludes the general framework. 



3. Results 

In what follows, we report results for a variety of different ensembles of sparse random 
matrices, in order to explore the capabilities and limitations of our approach. In order to 
properly appreciate the results presented below, it is worth pointing out that within our 
stochastic population-dynamics based approach to solving the fixed point equations fl30|) 
and (|3T|) . the integrals fl32l) . or (135]) . (1371) are evaluated by sampling from a population. 
Denoting by M the number of samples (a^, bi) taken, we have, e.g., 



as an approximation of (1371) . The e — > 0-limit is clearly singular in the first contribution to 
( 1391) . If bi + A 7^ for all bi in the sample, one obtains zero in the e — > 0-limit, whereas one 
obtains a diverging contribution, if bi + A = for at least one bi in the sample. The second 
alternative will quite generally be an event of probability zero, so a small regularizing 
e > must be kept in order to 'see' this contributions (if it exists). In what follows, we 
shall refer to the two contributions to (1371) . as p s (X) and p c (A), with 

a,j= a^>0 



Af 



Af 



A ) + -E 



7T 



ai + e 
{ ai + e) 2 + (bi + \y 



(39) 



a, : >0 
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The population-dynamics algorithm itself is run with a small regularizing e > (as 
required in (j2J) to guarantee existence of the integral). While running the algorithm, we 
use £ = 10~ 300 , which is close to the smallest representable real number in double-precision 
arithmetic on the machines used for the numerics. 



3.1. Poisson Random Graphs — Gaussian Couplings 

Our first results pertain to sparse matrices defined on Poisson random graphs, with 
Gaussian couplings. The left panel of Fig. [T] shows spectral densities for the case of mean 
connectivity c = 4, having Gaussian random couplings with (Kfj) = 1/c. For this system 
we find an integrable power-law divergence of the form 

p(A) ~ 0.05|A|-° 61 , A^O, (41) 

and a 5 peak at A = 0, the latter originating from isolated sites in the ensemble. Results 
of numerical diagonalizations (using a sample of 500 N x N matrices with N = 2000 are 
shown for comparison, and the agreement is excellent. 





Figure 1. Spectral density for matrices defined on Poissonian random graphs with c = 4 
(left panel) and c = 2 (right panel), having Gaussian random couplings with (Kfj) = 1/c. 
Full line: results obtained from the present theory; dashed line: results obtained from a 
sample of 2000x2000 matrices. In both cases e = 10 -300 was used in the evaluation of 

d&a. 



The behaviour changes rather drastically if the average connectivity is reduced to c = 2 
- a value closer to the percolation threshold c c = 1. In this case the spectral density shows 
strong fluctuations, when evaluated with the same small regularizer. These originate from 
pi in (BSD , and are related to the pure point spectrum associated with localized eigenstates 
coming from a collection of isolated finite clusters of all sizes in the ensemble. These 
exist for c = 4 as well, but their contribution is too small to be easily notable when 
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combined with ~p~ c in (I59"j) . In addition, there is a central 6 peak as in the c = 4-case which 
appears to be separated from the main bands by a gap; see the second panel in Fig [21 
The agreement with results of numerical diagonalization is fairly poor as it stands; in 
particular, exponential tails of localized states extending beyond the apparent edge of the 
central band are missed in this way. However, when (139]) is evaluated with a regularizing 
e = 10~ 3 comparable to the resolution of the A-scan, the agreement is once more excellent 
as shown in Fig [2j It is worth noting in this context that numerical simulations, in which 
binning of eigenvalues is used to determine the the spectral density also imply a form of 
regularization, and they do not distinguish continuous and singular contributions to the 
DOS if the distribution of the singular contributions is itself reasonably uniform. 

When displayed on a logarithmic scale, the results clearly reveal two interesting features: 
(i) a localization transition at A c ~ 2.295, characterised by a vanishing continuous 
contribution ~p~ c to (|39|) for |A] > A c , and (ii) exponential (Lifshitz) tails [19] in the spectral 
density, related to localized states represented by the singular contribution ~p~ s to (139])). and 
exhibited only through regularization. We shall substantiate this analysis in the following 
sub-section by looking at the behaviour inverse participation ratios. The same phenomena 
are seen for c = 4, where A c ~ 2.581. 

3.2. Inverse Participation Ratios and Localization 

In order to substantiate our identification of singular and continuous contributions to the 
spectral densities we look at Inverse Participation Ratios (IPRs) of eigenstates as obtained 
from numerical diagonalizations. Given eigenvectors v of a (random) matrix, their IPRs 
are defined as 



As eigenvectors can always be chosen to be normalized, we see that IPRs remain of order 
1 for localized states which have a few 0(1) eigenvector components — the extreme case 
being IPR(i?) = 1 for Vi = 5i^ — whereas they are Oi^N' 1 ) for fully extended states for 
which Vi = 0{N~ 1 / 2 ) for all i. 

Here we only produce a qualitative comparison for the two cases studied in the previous 
subsection, comparing IPRs computed for systems of size N = 100 and iV = 1000, and 
using scatter-plots of IPRs vs eigenvalues to exhibit the salient features. As clearly visible, 
there remains a substantial fraction of states at all A in the c = 2 case, which do not exhibit 
the iV -1 scaling of IPRs expected for delocalized states; the tails, and a small central band 
in particular appear to be dominated by localized states. By contrast in the c = 4 case 
there is a notable depletion of states with 0(1) IPRs, except for A = and in the tails of 




(42) 
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Figure 2. Upper left panel: Spectral density for matrices defined on Poissonian random 
graphs with c = 2 as in the previous figure, but now evaluated with a regularizing 
e = 1CP 3 in ([39| (full line). At the resolution given the result is indistinguishable from 
the numerical simulation results (dashed line) . Upper right panel: zoom into the central 
region comparing results obtained with the small regularizer, exhibiting a gap around 
the central peak (full line), with a larger regularizer e = 1CU 3 (short dashed line) and 
with results of numerical diagonalization (long dashed line). The same comparison is 
made in the lower panel for a larger portion of the spectrum on a logarithmic scale. 
The regularized e — 10~ 3 -results are on this scale indistinguishable from those of the 
numerical simulations. Note the localization transition and the Lifshitz tails as discussed 
in the main text. 

the spectrum. These findings are entirely consistent with our identifications made in the 
previous subsection. We note that the role of regularization in identifying localized states 
has been pointed out before using heuristics related to the evaluation of local densities of 
state [22] . 

We shall return to this issue in greater quantitative detail in a separate paper devoted to 
Anderson localization in discrete random Schrodinger operators defined on sparse random 
graphs [33J. 
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Figure 3. Scatterplots showing eigenvalue against IPRs for Poissonian random graphs 
with c = 2 (first row) and c = 4 (second row) . The graphs in the left column correspond 
to N = 100, those in the right column to N = 1000. 



3.3. Poisson Random Graphs — Bimodal Couplings 

We can also look at coupling distributions different from Gaussian for the non-zero 
couplings, e.g. fixed Kij = l/y/c or bi-modal = ±l/-y/c. As noted before [H], both 
give rise to the same spectral densities on large sparse (tree-like) graphs due to the absence 
of frustrated loops. It can also be seen as a consequence of the appearance of K 2 in (1251) . 

We choose a Poissonian random graph at the percolation threshold c = 1 as an example 
that allows us to highlight both the strengths and the limitations of the present approach. 
It is known that all states will be localized for this system. In Fig H] we compare results of 
a A-scan with resolution 5X = 10~ 3 , using a regularizer e = 10~ 4 for the scan. The smaller 
panels exhibit numerical diagonalization results, as well as a comparison between the two 
using a zoom into the region around A = 1. 

On the side of the strengths, we note that the spectral density obtained from our algorithm 
is able to display more details than can be exposed by simulation results obtainable at 
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100 : 




Figure 4. Comparison of spectral density for — ±l/y/c, on a Poissonian random 
graph with c = 1 as computed via the present algorithm (main panel) with results from 
numerical diagonalisation of JV x JV matrices of the same type with N — 2000 (lower 
left) and a direct comparison in the region around A = 1. 

reasonable effort. On the downside, one might note that the results for this system attain 
the status of semi-quantitative results, as they do depend on the chosen regularization, 
though in fairness it should be said that the same applies to the results obtained via 
numerical diagonalization where results vary with the binning resolution. In the present 
case this is due to the fact that the spectrum for most parts consists of a dense collection 
of 5 peaks [39]. A notable deficiency is the broadening of delta-peaks into Lorentzians of 
finite width, which creates artefacts around isolated delta-peaks, exemplified here by the 
peak at A = 0. Since the origin of this deficiency is understood, more precise details can, 
if desired, be recovered by choosing a smaller regularizing e. 
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3. 4- Regular and Scale- Free Random Graphs 

In the present section we consider matrices defined on regular and scale-free random 
graphs. 

3.4.1. Regular Random Graphs Our theory applies unmodified to matrices defined on 
graphs with degree distributions other than Poissonian, as long as the mean connectivity 
remains finite. We use this fact to obtain spectra of matrices with Gaussian random 
couplings defined on regular random graphs with fixed connectivity c, choosing (Kfj) = 
1/c for the couplings. Results for c = 4 and c = 100 are shown in Fig. The c = 4 results 
are compared with simulations, with results analogous to previous cases, including the 
presence of a localization transition at A c ~ 2.14 

The second example is chosen as a test to see the semicircular law jlQ] reemerge in the 
limit of large (though finite) connectivity. This limit can also be extracted from the fixed 
point equations. It is somewhat easier to verify for results pertaining to single instances 
[32J than for the ensemble. 




Figure 5. Spectral densities for a random graph with fixed connectivity c = 4 (left), 
and on a random graph with fixed non-random connectivity c = 100 (right). 



3.4.2. Scale-Free Graphs We have also looked at a scale free graph with connectivity 
distribution given by p(k) = -Po^ -7 with 7 = 4 and a lower cut-off at k — 2. Results 
shown in Fig. [6] reveal a continuous central band, and localized states for |A| > A c ~ 2.85 
much as in the other cases. For the present system, the tails in the spectral density follow 
a power law of the form p(A) ~ A 1-27 [Tfl |4"T] . 

Comparison with exact diagonalization results is facilitated by a fast algorithm that 
allows to generate sparse graphs with arbitrary degree distribution [42J. 
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Figure 6. Spectral density for for Kij = ±1/^ on a random graph with powerdaw 
degree distribution of average connectivity c ~ 2.623. Left panel: results obtained with 
small regularizer (full line), and numerical diagonalization results from a sample of 500 
matrices of dimension N — 2000 (dashed line). Right panel: the same results displayed 
on a logartithmic scale, this time with results regularized at e = 10~ 3 (short dashed line) 
included. 



3.5. Graph Laplacians 

Let us finally look at matrices row-constraints, such as related to discrete graph- 
Laplacians. 

The discrete graph Laplacian of a graph with connectivity matrix C = {%} has matrix 
elements 



Cij 5ij ^ ' c-ik . 



(43) 



A quadratic form involving the Laplacian can be written in the form 

1 2 ^ijUiUj = — CijiUi - Ujf . 



2 t-^i -^ij U i U j i 
ij 



(44) 



As before we shall be interested in more general matrices with zero row-sum constraint 
of the form 



■Mij Cij Kij 5ij ^ ' Ci^Kij 
k 



(45) 



To evaluate the spectral density within the present framework one would thus have to 
compute 



K 
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instead of The required modification has, of course, been noted earlier [T51 |4"3"] . 
The resulting problem constitutes precisely (the harmonic variant of) the translationally 
invariant systems, for which the framework in [30] was developed in the first place. The 
general theory can be copied word for word, and the fixed point equations ( 1301) . ( I3TT) 
remain formally unaltered except for the change in Z2(u;,u',K) in (|23|) . owing to the 
modified interaction term, which gives rise to a modified expression for Cl(uj', K) in (1281) . 
We obtain 

Ki / 

Ci(u/, K) = (46) 

instead of (1281 . Fig. [7] shows the spectrum of a Laplacian for a Posisson random graph with 
c = 2, comparing our solution (upper left panel) computed with e = 10 -3 with numerical 
diagonalization results in the upper right panel. We use = 1/c for the non-zero matrix 
elements in this case. As in the other cases, we observe a localization transition, here at 
A c ~ —3.98. Results obtained with a small regularizer e = 10~ 300 exhibiting only the 
continuous part of the spectrum are shown in the lower panel. 



3.6. Unfolding Spectral Densities 

As a last item in this study we look at the possibility of unfolding the spectral density 
according to contributions of local densities of state, coming from vertices of different 
coordination, as suggested by Eq. ( |32l) . This method has been used in [30] to look at 
distributions of Debye- Waller factors in amorphous systems, unfolded according to local 
coordinations. In the present context it may provide an interesting diagnostic tool to help 
understanding localization phenomena. 

Fig [8] exhibits the spectrum of the graph Laplacian shown in the previous figure along 
with its unfolding into contributions of local densities of state with different coordination. 
The present example clearly shows that — somewhat paradoxically — the well connected 
sites are the ones providing the dominant contributions to localized states in the lower 
band-edge Lifshitz tails. The clearly identifiable humps in the figure correspond from left 
to right to k = 9, k = 8, k = 7, k = 6, k = 5, k = 4, and k = 3, which easily allows 
to identify the corresponding contributions to the spectral density, the contribution of 
k = 2 gives rise to several notable humps in the spectral density, and together with the 
k = 1 contribution is mainly responsible for the dip at A = —1. The k = contribution 
is mainly responsible for the 5-peak at A = (which is broadenend into a Lorentzian of 
width e = 10 -3 due to the regularization, as discussed earlier. 
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Figure 7. Spectral density for the Laplacian on a Poissonian random graph with c = 2 
as computed via the present algorithm. Upper left panel: e = 10~ 3 -results; upper right 
panel: results from numerical diagonalisation of N x N matrices of the same type with 
N = 2000. Lower panel: continuous part of the spectrum obtained using e = 10~ 300 as a 
regularizer. 

4. Conclusions 

In the present paper we have used a reformulation of the replica approach to 
the computation of spectral densities for sparse matrices, which allows to obtain 
spectral densities in the thermodynamic limit to any desired detail — limited only by 
computational resources. Our method is versatile in that it allows to study systems with 
arbitrary degree distributions, as long as they give rise to connectivity distributions with 
finite mean. A cavity approach that emphasises results on finite instances will appear 
elsewhere [32]. As expected (and well known), the Wigner semi-circle reemerges in the 
large c limit as discussed in |32j . Large and small A asymptotics remain to be investigated. 
Our method allows to expose the separate contributions of localized and extended states 
to the spectral density, and thereby to study localization transitions. We shall explore this 
issue in greater detail in a separate publication. Indeed, with results for graph-Laplacians 
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Figure 8. Spectral density for the Laplacian on a Poissonian random graph with c = 2 
(full upper line) , shown together with its unfolding according to contributions of different 
coordination, as discussed in the main text. 

in hand, the step towards a study of discrete random Schrodinger operators and Anderson 
localization in such systems is just around the corner [33]. A generalization to asymmetric 
matrices using both the cavity method and a replica approach for the ensemble along the 
lines of [H] is currently under investigation in our group [33]. Other problems we have 
started to look at are spectra of modular systems [46] and small world networks. 

We believe our results to constitute an improvement over previous asymptotic results as 
well as over results obtained by closed form approximations. They may open the way to 
further interesting lines of research. Let us here mention just a few such examples: within 
RMT proper, one might wish to further investigate the degree of universality of level 
correlations in these systems [UJ ; one could refine the random matrix analysis of financial 
cross-correlations [7J by taking non-trivial degree distributions of economic interactions 
into account, or one might wish to look at finite connectivity variants of random reactance 
networks [IE], taking e.g. regular connectivity 4 to compare with results of numerical 
simulations of such systems on two-dimensional square lattices. 

Acknowledgements It is a pleasure to thank Giiler Ergiin, Jort van Mourik, Isaac Perez- 
Castillo, Tim Rogers and Koujin Takeda for illuminating discussions. Jort van Mourik also 
kindly supplied instances of scale free-random graphs to allow comparison of ensemble 
results and results from numerical diagonalization in this case. 
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Appendix A. Population Dynamics 

The stochastic algorithm used to solve (130|) . (13T1) takes the following form. Populations 
{^f, 1 < i < A/p} and {ujf, 1 < z < AT p } are randomly initialized with Re cUj > and 
Re Cji > 0. 

Then the following steps are iterated 

1. Generate a random k ~ \p c {k)- 

2. Randomly select — 1 elements from {cc^; 1 < i < N p }; compute 

fe-i 

fi = iA e + X)wi J , (A.l) 
i=i 

and replace by f2 for a randomly selected i G {1, . . . , N p }. 

3. Select j G {1, . . . , Ap} at random, generate a random K according to distribution of 

bond strengths; compute 

Cl = , (or (l = — for zero row-sums ] , (A.2) 

and replace Cji by Cl for a randomly selected i G {l,...,Ap}. 

4. return to 1. 

This algorithm is iterated until populations with stable distributions of {Cji; 1 < i < N p } 
and {u>i] 1 < i < N p } are attained. 

A variant of this algorithm when implemented on instances of real graphs generates the 
belief-propagation or cavity equations for this problem, as studied in [32]. It can be derived 
directly in terms iterative evaluations of ([2]) on locally tree-like graphs. 
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